!! Copyright (C) 2009,2010,2011,2012  Marco Restelli
!!
!! This file is part of:
!!   LDGH -- Local Hybridizable Discontinuous Galerkin toolkit
!!
!! LDGH is free software: you can redistribute it and/or modify it
!! under the terms of the GNU General Public License as published by
!! the Free Software Foundation, either version 3 of the License, or
!! (at your option) any later version.
!!
!! LDGH is distributed in the hope that it will be useful, but WITHOUT
!! ANY WARRANTY; without even the implied warranty of MERCHANTABILITY
!! or FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public
!! License for more details.
!!
!! You should have received a copy of the GNU General Public License
!! along with LDGH. If not, see <http://www.gnu.org/licenses/>.
!!
!! author: Marco Restelli                   <marco.restelli@gmail.com>


!>\brief
!!
!! Nonlinear, breaking mountain wave.
!!
!! \n
!!
!! For this test case we refer to Chapt. 5.2.6 of [Restelli M., PhD.
!! thesis, Milan 2007], and in particular to the breaking wave case of
!! Fig. 5.21.
!<----------------------------------------------------------------------
module mod_breaking_wave_test

!-----------------------------------------------------------------------

 use mod_messages, only: &
   mod_messages_initialized, &
   error,   &
   warning, &
   info

 use mod_kinds, only: &
   mod_kinds_initialized, &
   wp

 use mod_fu_manager, only: &
   mod_fu_manager_initialized, &
   new_file_unit

 use mod_physical_constants, only: &
   mod_physical_constants_initialized, &
   t_phc

!-----------------------------------------------------------------------
 
 implicit none

!-----------------------------------------------------------------------

! Module interface

 public :: &
   mod_breaking_wave_test_constructor, &
   mod_breaking_wave_test_destructor,  &
   mod_breaking_wave_test_initialized, &
   test_name,        &
   test_description, &
   phc,              &
   ntrcs,            &
   ref_prof, t_ref,  &
   coeff_dir,        &
   coeff_neu,        &
   coeff_norm,       &
   coeff_sponge,     &
   coeff_visc,       &
   coeff_init,       &
   coeff_outref

 private

!-----------------------------------------------------------------------

 ! public interface
 character(len=*), parameter   :: test_name = "breaking-wave"
 character(len=100), protected :: test_description(5)
 type(t_phc), save, protected  :: phc
 integer, parameter            :: ntrcs = 0
 character(len=*), parameter   :: ref_prof = "isothermal"
 real(wp), target              :: t_ref

 logical, protected :: mod_breaking_wave_test_initialized = .false.
 ! private members
 real(wp) :: &
   nu,     & ! kinematic viscosity
   nn,     & ! Brunt-Vaisala freq.
   t_s,    & ! surface temperature
   u_mean, & ! mean flow
   ! mountain geometry
   hm, a, x0, &
   ! sponge layers
   zd, zt, alpha, dv, dh, x1(2), x2(2)
 character(len=*), parameter :: &
   test_input_file_name = 'breaking_wave.in'
 character(len=*), parameter :: &
   this_mod_name = 'mod_breaking_wave_test'

 ! Input namelist
 namelist /input/ &
   nu, nn, t_s, u_mean, hm, a, x0, zd, zt, alpha, dv, dh, x1, x2

!-----------------------------------------------------------------------

contains

!-----------------------------------------------------------------------

 subroutine mod_breaking_wave_test_constructor(d)
  integer, intent(in) :: d

  integer :: fu, ierr
  character(len=*), parameter :: &
    this_sub_name = 'constructor'

   !Consistency checks ---------------------------
   if( (mod_messages_initialized.eqv..false.)   .or. &
       (mod_kinds_initialized.eqv..false.)      .or. &
       (mod_fu_manager_initialized.eqv..false.) .or. &
       (mod_physical_constants_initialized.eqv..false.) ) then
     call error(this_sub_name,this_mod_name, &
                'Not all the required modules are initialized.')
   endif
   if(mod_breaking_wave_test_initialized.eqv..true.) then
     call warning(this_sub_name,this_mod_name, &
                  'Module is already initialized.')
   endif
   !----------------------------------------------

   ! Read input file
   call new_file_unit(fu,ierr)
   open(fu,file=trim(test_input_file_name), &
      status='old',action='read',form='formatted',iostat=ierr)
    if(ierr.ne.0) call error(this_sub_name,this_mod_name, &
      'Problems opening the input file')
    read(fu,input)
   close(fu,iostat=ierr)

   write(test_description(1),'(a,i1)') &
     'Breaking wave, dimension d = ',d
   write(test_description(2),'(a,e9.3,a)') &
     '  viscosity nu = ',nu,';'
   write(test_description(3),'(a,e9.3,a)') &
     '  mean flow velocity = ',u_mean,';'
   write(test_description(4),'(a,e9.3,a)') &
     '  Brunt-Vaisala frequency = ',nn,';'
   write(test_description(5),'(a,e9.3,a)') &
     '  surface temperature = ',t_s,'.'

   ! All the physical constants use the default values: phc unchanged

   t_ref = 300.0_wp

   mod_breaking_wave_test_initialized = .true.
 end subroutine mod_breaking_wave_test_constructor

!-----------------------------------------------------------------------
 
 subroutine mod_breaking_wave_test_destructor()
  character(len=*), parameter :: &
    this_sub_name = 'destructor'
   
   !Consistency checks ---------------------------
   if(mod_breaking_wave_test_initialized.eqv..false.) then
     call error(this_sub_name,this_mod_name, &
                'This module is not initialized.')
   endif
   !----------------------------------------------

   mod_breaking_wave_test_initialized = .false.
 end subroutine mod_breaking_wave_test_destructor

!-----------------------------------------------------------------------

  !> Viscosity
  pure function coeff_visc(x,uu) result(mu)
   real(wp), intent(in) :: x(:,:), uu(:,:)
   real(wp) :: mu(size(x,2),1+size(x,1))

    mu = nu
  end function coeff_visc

!-----------------------------------------------------------------------

 !> Constant stability atmosphere with uniform horizontal flow.
 pure function coeff_init(x,u_r) result(u0)
  real(wp), intent(in) :: x(:,:), u_r(:,:)
  real(wp) :: u0(2+size(x,1),size(x,2))

  real(wp) :: delta
  real(wp), dimension(size(x,2)) :: z, zexp, theta, pi, p, t

   ! first define a stable, constant stability atmosphere
   z = x(size(x,1),:) ! vertical coordinate
   delta = 1.0_wp - phc%gravity**2/(phc%cp*nn**2*t_s)
   zexp = exp( nn**2/phc%gravity * z )
   theta = t_s * zexp
   pi = (1.0_wp-delta)/zexp + delta

   ! translate into primitive variables
   p = phc%p_s*pi**(1.0_wp/phc%kappa)
   t = pi*theta

   ! finally translate the result in conservative variables
   u0(1,:) = p/(phc%rgas*t)                                        ! rho
   u0(2,:) = u0(1,:)*(phc%cv*t + 0.5_wp*u_mean**2 + phc%gravity*z) ! e
   u0(3,:) = u0(1,:)*u_mean                                        ! Ux
   u0(4:,:) = 0.0_wp                                               ! Uyz
   ! subtract the reference state
   u0(1,:) = u0(1,:) - u_r(1,:) ! rho
   u0(2,:) = u0(2,:) - u_r(2,:) ! e

 end function coeff_init
 
!-----------------------------------------------------------------------

 !> Use the initial condition.
 pure function coeff_outref(x,u_r) result(u0)
  real(wp), intent(in) :: x(:,:), u_r(:,:)
  real(wp) :: u0(2+size(x,1),size(x,2))

   u0 = coeff_init(x,u_r)

 end function coeff_outref

!-----------------------------------------------------------------------

 !> Use the initial condition.
 pure function coeff_dir(x,u_r,breg) result(ub)
  real(wp), intent(in) :: x(:,:)
  real(wp), intent(in), optional :: u_r(:,:)
  integer, intent(in), optional :: breg
  real(wp) :: ub(2+size(x,1),size(x,2))

   ub = coeff_init(x,u_r)

 end function coeff_dir

!-----------------------------------------------------------------------

 pure function coeff_neu(x,u_r,breg) result(h)
  real(wp), intent(in) :: x(:,:)
  real(wp), intent(in), optional :: u_r(:,:)
  integer, intent(in), optional :: breg
  real(wp) :: h(size(x,1),2+size(x,1),size(x,2))

   h = 0.0_wp

 end function coeff_neu

!-----------------------------------------------------------------------

 !> Normal to the boundary
 !!
 !! \warning This function should be only called for the bottom
 !! boundary.
 !!
 !! In 2D, we have
 !! \f{equation}{
 !!  y(x) = \frac{h_m}{1+\left(\frac{x-x_0}{a}\right)^2}
 !! \f}
 !! so that
 !! \f{equation}{
 !!  \tilde{\underline{n}}(x) = 
 !!   \left[\begin{array}{c} \frac{dy}{dx}(x) \\\ -1 \end{array} \right]
 !!   = \left[\begin{array}{c} 
 !!    -2\frac{h_m}{a}\frac{\xi}{\left(1+\xi^2\right)^2} \\\ -1
 !!   \end{array} \right]
 !! \f}
 !! where the \f$\tilde{}\f$ indicates that the resulting vector must be
 !! normalized and
 !! \f{equation}{
 !!  \xi = \frac{x-x_0}{a}.
 !! \f}
 !<
 pure function coeff_norm(x,breg) result(nb)
  real(wp), intent(in) :: x(:,:)
  integer, intent(in) :: breg
  real(wp) :: nb(size(x,1),size(x,2))

  integer :: l
  real(wp) :: xi(size(x,2))

   breg_case: select case(breg)
    case(3) ! bottom boundry, 2D
     xi = (x(1,:) - x0)/a
     nb(1,:) = -2.0_wp*hm/a * xi/(1.0_wp+xi**2)**2
     nb(2,:) = -1.0
     do l=1,size(nb,2)
       nb(:,l) = nb(:,l)/sqrt(nb(1,l)**2+nb(2,l)**2)
     enddo
    case default
     nb = huge(1.0_wp) ! this is an error
   end select breg_case

 end function coeff_norm

!-----------------------------------------------------------------------

 pure subroutine coeff_sponge(tmask,tau,taup,x)
  real(wp), intent(in) :: x(:,:)
  logical, intent(out) :: tmask
  real(wp), intent(out) :: tau(:), taup(:)

  real(wp), parameter :: pi_trig = 3.141592653589793_wp, &
    big = huge(1.0_wp)
  integer :: l, i
  real(wp) :: z, zn, d_bnd, coeff

   ! initialization
   tau = 0.0_wp
   taup = 0.0_wp
   tmask = .false.

   do l=1,size(x,2)

     ! upper layer
     z = x(size(x,1),l)
     if(z.gt.zd) then ! upper layer nonzero

       tmask = .true.

       ! gravity wave damping
       zn = (z-zd)/(zt-zd) ! normalized vertical coordinate
       if(zn.le.0.5_wp) then
         tau(l) = tau(l) + alpha/2.0_wp*(1.0_wp-cos(zn*pi_trig))
       else
         tau(l) = tau(l) + alpha/2.0_wp*(1.0_wp+(zn-0.5_wp)*pi_trig)
       endif

       ! acoustic wave damping
       d_bnd = (zt-z)/dv ! normalization
       if(d_bnd.le.3.0_wp) then
         ! tmask already set to .true.
         if(d_bnd.lt.1000.0_wp/big) then
           coeff = big/1000.0_wp
         else
           coeff = (1.0_wp-tanh(d_bnd))/tanh(d_bnd)
         endif
         tau(l)  = tau(l)  + coeff
         taup(l) = taup(l) + coeff
       endif

     endif

     ! lateral sponges (only acoustic wave damping)
     d_bnd = big
     do i=1,size(x,1)-1 ! first d-1 coordinates
       d_bnd = min( d_bnd , min( abs(x(i,l)-x1(i)) , abs(x(i,l)-x2(i)) ) )
     enddo
     d_bnd = d_bnd/dh ! normalization
     if(d_bnd.le.3.0_wp) then
       tmask = .true.
       if(d_bnd.lt.1000.0_wp/big) then
         coeff = big/1000.0_wp
       else
         coeff = (1.0_wp-tanh(d_bnd))/tanh(d_bnd)
       endif
       tau(l)  = tau(l)  + coeff
       taup(l) = taup(l) + coeff
     endif

   enddo
   
 end subroutine coeff_sponge

!-----------------------------------------------------------------------

end module mod_breaking_wave_test

